#include "least_square_fitting.h"
#include <cmath>
#include <cstdlib>

namespace rmd {

CLeastSquare::CLeastSquare() {}

CLeastSquare::~CLeastSquare() {}

/*************************************************************************
系数矩阵的限幅处理，防止它溢出，目前这个函数很不完善，并不能很好地解决这个问题
原理：矩阵解行列式，同一行乘以一个系数，行列式的解不变
以行为单位处理
**************************************************************************/
int CLeastSquare::ParalimitRow(double* para, int sizeSrc, int row) {
  double temp;
  double max = abs(ParaBuffer(para, row, 0));
  double min = max;
  for (int i = sizeSrc; i > 0; --i) {
    temp = abs(ParaBuffer(para, row, i));
    if (max < temp) {
      max = temp;
    }
    if (min > temp) {
      min = temp;
    }
  }

  double factor = (max + min) * 0.000005 + 0.0000001;
  for (int i = sizeSrc; i >= 0; --i) {
    ParaBuffer(para, row, i) /= factor;
  }

  return 0;
}

/****************************************************************************
同上，以矩阵为单位处理
****************************************************************************/
int CLeastSquare::Paralimit(double* para, int sizeSrc) {
  for (int i = 0; i < sizeSrc; ++i) {
    ParalimitRow(para, sizeSrc, i);
  }
  return 0;
}

/****************************************************************************
系数矩阵行列式变换
****************************************************************************/

int CLeastSquare::ParaPreDealA(double* para, int sizeSrc, int size) {
  size -= 1;
  for (int i = 0; i < size; ++i) {
    for (int j = 0; j < size; ++j) {
      ParaBuffer(para, i, j) =
          ParaBuffer(para, i, j) * ParaBuffer(para, size, size) -
          ParaBuffer(para, size, j) * ParaBuffer(para, i, size);
    }
    ParaBuffer(para, i, sizeSrc) =
        ParaBuffer(para, i, sizeSrc) * ParaBuffer(para, size, size) -
        ParaBuffer(para, size, sizeSrc) * ParaBuffer(para, i, size);
    ParaBuffer(para, i, size) = 0;
    ParalimitRow(para, sizeSrc, i);
  }

  return 0;
}

/*****************************************************************************
系数矩阵行列式变换，与ParaPreDealA配合
完成第一次变换，变换成三角矩阵
*****************************************************************************/
int CLeastSquare::ParaDealA(double* para, int sizeSrc) {
  for (int i = sizeSrc; i > 0; --i) {
    ParaPreDealA(para, sizeSrc, i);
  }

  return 0;
}

/*****************************************************************************
系数矩阵变换
*****************************************************************************/
int CLeastSquare::ParaDeal(double* para, int sizeSrc) {
  Paralimit(para, sizeSrc);
  ParaDealA(para, sizeSrc);
  ParaDealB(para, sizeSrc);

  return 0;
}

/******************************************************************************
系数矩阵行列式变换
******************************************************************************/
int CLeastSquare::ParaPreDealB(double* para, int sizeSrc, int OffSet) {
  for (int i = OffSet + 1; i < sizeSrc; ++i) {
    for (int j = OffSet + 1; j <= i; ++j) {
      ParaBuffer(para, i, j) *= ParaBuffer(para, OffSet, OffSet);
    }
    ParaBuffer(para, i, sizeSrc) =
        ParaBuffer(para, i, sizeSrc) * ParaBuffer(para, OffSet, OffSet) -
        ParaBuffer(para, i, OffSet) * ParaBuffer(para, OffSet, sizeSrc);
    ParaBuffer(para, i, OffSet) = 0;
    ParalimitRow(para, sizeSrc, i);
  }

  return 0;
}

/******************************************************************************
系数矩阵行列式变换，与ParaPreDealB配合
完成第一次变换，变换成对角矩阵，变换完毕
******************************************************************************/
int CLeastSquare::ParaDealB(double* para, int sizeSrc) {
  for (int i = 0; i < sizeSrc; ++i) {
    ParaPreDealB(para, sizeSrc, i);
  }

  for (int i = 0; i < sizeSrc; ++i) {
    if (ParaBuffer(para, i, i)) {
      ParaBuffer(para, i, sizeSrc) /= ParaBuffer(para, i, i);
      ParaBuffer(para, i, i) = 1.0;
    }
  }

  return 0;
}

/******************************************************************************
最小二乘法的第一步就是从XY数据里面获取系数矩阵
double* para：系数矩阵地址
const double* x：X数据地址
const double* y：Y数据地址
int amount：XY数据组数
int sizeSrc：系数矩阵大小（sizeSrc）行（sizeSrc + 1）列
******************************************************************************/

int CLeastSquare::GetParaBuffer(double* para, const double* x, const double* y,
                                int amount, int sizeSrc) {
  for (int i = 0; i < sizeSrc; ++i) {
    ParaBuffer(para, 0, i) = 0;
    int i_index = 2 * (sizeSrc - 1) - i;
    for (int j = 0; j < amount; ++j) {
      ParaBuffer(para, 0, i) += pow(*(x + j), i_index);
    }
  }

  for (int i = 1; i < sizeSrc; ++i) {
    ParaBuffer(para, i, sizeSrc - 1) = 0;
    int i_index = (sizeSrc - 1) - i;
    for (int j = 0; j < amount; ++j) {
      ParaBuffer(para, i, sizeSrc - 1) += pow(*(x + j), i_index);
    }
  }

  for (int i = 0; i < sizeSrc; ++i) {
    ParaBuffer(para, i, sizeSrc) = 0;
    int i_index = (sizeSrc - 1) - i;
    for (int j = 0; j < amount; ++j) {
      ParaBuffer(para, i, sizeSrc) += (*(y + j)) * pow(*(x + j), i_index);
    }
  }

  for (int i = 1; i < sizeSrc; ++i) {
    for (int j = 0; j < sizeSrc - 1; ++j) {
      ParaBuffer(para, i, j) = ParaBuffer(para, i - 1, j + 1);
    }
  }

  return 0;
}

/*****************************************************************************
整个计算过程
******************************************************************************/
int CLeastSquare::Cal(const double* bufferX, const double* bufferY, int amount,
                      int sizeSrc, double* paraResK) {
  double* paraK = (double*)malloc(sizeSrc * (sizeSrc + 1) * sizeof(double));
  GetParaBuffer(paraK, bufferX, bufferY, amount, sizeSrc);
  ParaDeal(paraK, sizeSrc);
  for (amount = 0; amount < sizeSrc; ++amount) {
    *paraResK = ParaBuffer(paraK, amount, sizeSrc);
    ++paraResK;
  }

  free(paraK);
  return 0;
}
}
